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Abstract 

Molecular process of polymer collapse was reproduced by isothermal molec- 
ular dynamics simulation. The initial polymer chains were obtained by mean 
of random walks in continuum space. Two potential models were considered to 
represent short range interactions between monomers. Various structural prop- 
erties during the collapse of the polymer were measured and different collapse 
pathways were observed. For low temperatures, we have obtained a spontaneous 
collapse to helical states. 
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1. Introduction 

During the last four decades, the collapse of polymer chains has been exten- 
sively investigated [1~3]. The key reason for the interest to polymer collapse is 
the close relationship of this process to the protein folding, one of the most chal- 
lenging problems in molecular biology. It is well known that in a good solvent, 
a polymer chain is extended. However, as one decrease the solvent quality, the 
polymer chain adopts a conformation which tries to minimize solvent-polymer 
contacts and the final shape of the chain is compact [4]. The collapse of poly- 
mers has been studied by many authors recently using both analytical methods 
and computer simulations [5-7] . They all strongly suggest that the polymer col- 
lapse is a multistage process. Schnurr, Mackintosh and Williams [8-10] treated 
the case of a single semiflexible polymer in solution by means of Brownian Dy- 
namics. They found that the semiflexible polymer collapses to a torus and the 
globular ground state is unfavorable. Frish and Verga [11, 12] studied the ef- 
fect of the solvent on the collapse pathways of a single flexible homopolymer 
chain. The results found by this group show that the collapse of an isolated 
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homopolymcr is dominated by the second order coil to globule transition. De- 
pending on the quench dept they observed intermediate sausage regime or a set 
of pearls along the chain followed by shrink in polymer size to a final spherical 
globule. These results agree with the general scenario for the kinetics of the 
collapse of a flexible polymer coil presented by de Gennes [4]. Although, fully 
explicit simulations with solvent particles are still very expensive numerically, 
mainly in three dimensions. Poison et al [13 15] have made in a series of works 
a molecular dynamics simulation of a heteropolymer immersed in an explicitly 
modeled solvent. They studied the effect of hydrophobic, hydrophilic monomers 
distribution on the collapse process. A more recent study [16], investigated the 
role of hydrodynamics on the kinetics of polymer collapse transition. This study 
shows that hydrodynamics interactions speed up the collapse of the polymer and 
enhance cooperative motion of the monomers. 

The purpose of the present work is to study the effect of the temperature 
and the short range interactions potential model on the collapse pathways of 
polymers using molecular dynamics simulation. 

Simulating protein folding remains challenging and the question is how can 
one reduce computational efforts and correctly predict the folding kinetics? 

The paper is organized as follows. In section (2) we describe the details of 
the simulation method used in the present study. In section (3) we report our 
results for the structural properties during the process of the polymer collapse 
than we conclude with a discussion of the results. 

2. Theory and methods 

2.1 Molecular model 

We have used the bead spring model to represent a polymer chain of N 
monomers. Non bonded monomers interact via a truncated and shifted Lennard- 
Jones potential 



Where C (re) ensures that the potential vanishes at the cut-off distance 
Tc = 2.5. We have used standard, reduced units in our calculations, where 
all distances and energies are expressed in terms of Lennard-Jones parameters 
a and £, respectively; mass is expressed in terms of the monomer mass m and 
the unit of time is ^J ma^ je. Temperature and energy are made numerically 
identical by setting the Bolzmann constant fcs to unity. 

We consider two different potential models for representing the nearest neigh- 
bor monomers connectivity, henceforth referred to as potential models A and B. 
In potential model A, monomer i interacts with monomer i — 1 and monomer 
i + 1 via a strong anharmonic potential defined by 





else 



(1) 



\JA{r) = a{r - r^f 4- fo (r - tq)'^ 



(2) 
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Here To = 0.85 is the equilibrium distance. We take a = 30 and b = 100 
for the constants characterizing the harmonic and anharmonic interactions as 
in [11,12]. In case of the potential model B, the interactions between pairs of 
bonded monomers are represented only by a harmonic potential of the form 

UB{r) = a{r-rof (3) 

2.2 Initial Configuration 

The initial configuration of the polymer chain is a three dimensional random 
walk in continuum space . Starting from a position at the origin ri (0, 0, 0), the 
positions {xi, yi,Zi) of the monomers are generated using the following procedure 

Xi = Xi-i + acos{<p)sin{6) 

< yi = yi-i + asin{(j>)sin{9) (4) 

Zi = Zi_i+ acos{6) 

9 and (f) are two uniformly random angles chosen respectively in the ranges 
[0, 7r/9] and [0, 27r] and a is the distance between two successive monomers. 9 
has to be small to simulate a stretched polymer chain. In Figure (1), we show 
a snapshot of a polymer chain generated by this procedure. 

2.3 Polymer Dynamics 

We have performed molecular dynamics simulations in the canonical ensem- 
ble. In order to maintain the temperature constant, we have coupled the system 
to the Nose Hoover Thermostat [17, 18]. 

This scheme derives from a modified Lagrangian formalism that introduces 
a new degree of freedom. The equations of motion are then 

r'i = Pi /mi 

< Pi = -^^-^Pi (5) 
. ^ = ^[Ef=iS-3iVfcBT] 

Here and Pi are respectively the position and momentum of monomer i. Q 
is a parameter defining thermal inertia in the system, it must be carefully chosen 
to prevent the simulation from unwanted oscillations. In our simulations, we 
have chosen Q = 10. Kb represents the Boltzman constant and ^ is the friction 
variable. We have solved numerically the equations of motion, using the Verlet 
Newton-Raphson algorithm [19] and the time step was set to At = 10~^. 
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2-4 Ohservahles 

The quantities of central interest sampled during the simulations are the 
gyration radius Rg of the polymer and the long range order parameter S. For 
a single chain conformation {tn}, they are defined respectively by 



Rg - 



i 



1 ^ 



(6) 



where is the position of the ith monomer, N is the number of monomers, 
and rem. is the center of mass position of the polymer chain given by Tc.m. = 
W E»=i n and 



^ JV-2 



(7) 



The parameter S measures the long-range order present in the folded struc- 
tures. 



2.5 Locating the @ -point 

From theorical point of view, the 0-point is the temperature at which the 
distribution of the end-to-end distance of the polymer chain P{R) display gaus- 
sian statistics 

P{R) = AnR'^Aiexpi-TrAR'^) (8) 
Where A = 3/(27r(ii^)) and R is the end-to-end vector defined by 

N 

R = Y,bi (9) 

i=l 

Here, bi is the bond vector. 

We have used a method described by Yong et al [20] and based on Rosen- 
bluth sampling [21] for locating the 6-point. The main idea of this method 
is to compute the distribution of the end to end distance and compare it to 
the theorical distribution P{R) for different values of the temperature. In our 
simulations, we have explored the temperature range T — 3.0 — 4.0 and found 
the 8-point at the temperature ~ 3.4 for the set of parameters used in the 
molecular model presented above. We show in figure (2) that for a polymer 
chain size N = 256, the end-to-end distance distribution fit closely the theorical 
distribution for the temperature w 3.4. 
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3. Results and discussion 

Wc have investigated the collapse pathways of the polymer chains after an 
abrupt decrease in temperature. We first consider the case of small quench 
below the 8-point T = 3.0. Immediately after the quench, pearls start to 
grow along the polymer chain backbone, these pearls nucleate more easily from 
the polymer chain extremities as described in [11]. The high mobility of the 
monomers helps the pearls to merge rapidly to a compact globule. Figures (3) 
and (4) show a representative picture of this process. It is important to note 
that for a small temperature quench, similar collapse pathways were observed 
during the simulations cither using the potential models A or B. 

In the case of a strong quench T = 0.1, the polymer collapse pathway is 
completely different. We have observed a well ordered helix structure along the 
chain. Figure (5) and (6) show the spontaneous collapse of the polymer chain 
into this structure. (The polymer is drawn as a tube to see clearly the helix 
structure) . 

The rate of collapse kinetics was measured by monitoring the variation of the 
gyration radius Rg with time. Figure (7) shows that the deeper the temperature 
quench, the longer the collapse. 

We found that in simulations using potential model A, polymer chains show 
a better folding behaviors and a more pronounced helical ground states than 
in simulations using potential model B, in agreement with [22], where it was 
argued that a plain harmonic potential could induce energy localization in some 
specific modes and significantly increase the time for equilibration. This result 
is clearly visible in Figure (8) through the plot of the long range order parameter 
5 as a function of time for both potential models A and B. 

The numerical simulation of the collapse pathways of polymer chains in the 
present work and its comparison with previous studies can be summarized as 
follows: 

(1) The process of collapse of polymer chains is a two stage process, in agree- 
ment with the results found by Frish and Verga [11]. For a small temperature 
quench below the ©-point, the first stage is characterized by the formation of 
pearls along the chain and in the second stage, the polymer becomes compact. 
For a strong temperature quench, the polymer collapses spontaneously to an 
ordered helix structure without the use of any torsional potential as in [23, 24] 
or a complicated potential form as in [25] . Special ground states like double he- 
lices with antiparallel alignment have been also obtained at some temperature 
quenches (see figure (9)). 

(2) The polymer chain collapse seems to be much slower for a strong tem- 
perature quench below the 0-point. 

(3) The use of potential model A shows a better ordered structures at low 
temperatures. 

The simplified approach of this work has been used in a polymer context, 
however it is also relevant to protein folding as the helical states are the main 
ground states observed experimentally [26]. A better investigation of those 
ground states would be a future possible direction. The role of monomers rel- 
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ative contact order on the spontaneous polymer collapse into helical states is 
under investigation. 
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Figure 1: Initial configuration of the polymer chain 
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Theorical end-to-end distance distribution 

■ Caicuiated end-to-end distance distribution 




Figure 2: End-to-end distance distribution for the temperature 9 = 3.4 




Figure 3: Molecular dynamics simulation of a iV = 256 polymer chain at tem- 
perature T = 3.0, pearls regime at f = 4 x 10^ At 
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Figure 4: Molecular dynamics simulation of a N = 256 polymer chain at tem- 
perature T = 3.0, globular state a,t t = 10^ At 
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Figure 6: Molecular dynamics simulation of a A?^ = 256 polymer chain at tem- 
perature T = 0.1, a well formed helix at t = 3 x 10^ At 
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Figure 7: Gyration radius Rg versus time for two differents temperatures T 
0.1 and T = 3.0 and polymer chain size A'' = 256 
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— ° — Potential model A 

— o — Potential model B 
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Figure 8: Order parameter as function of time for polymer chain size N = 256 
at temperature T = 0.1 using potential model A and potential model B 
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